*Main Estimates for "Reveal Preference from Voluntary Contribution" (Bobolink SBC)

clear all
set more off
capture log close

*Specify directories here
global dta "C:\Users\zhenshanchen\Documents\RP from VC\Valuing public good WTP\dta"
global results "C:\Users\zhenshanchen\Documents\RP from VC\Valuing public good WTP\Results"

set seed 123456789

 
**********************************
* Summary Statistics -  Table 3  *
**********************************
clear all
set more off
use "$dta\reanalysis_hurdle_final0.dta",clear
replace MinAmt=. if dc==1
drop if dcamt>offeramt&offeramt_all>0&dc==1
global X "donate acres view numboboter NFledglings_recode MinAmt dcamt dc PM UPA UPC round2"
eststo BoboSS_dc: estpost sum $X if dc==1
eststo BoboSS_oe: estpost sum $X if dc==0
esttab BoboSS_dc BoboSS_oe using "$results\SummaryStatistics1.csv",replace cells("mean(fmt(3))" "sd(fmt(3) par)")

global HH "ppower age female dany denv moany momags moelect mokids round2 Refund07 southto138 lnparceldens PObox lncoastdist"
eststo BoboSS_hh_dc: estpost sum $HH if dc==1
eststo BoboSS_hh_oe: estpost sum $HH if dc==0
esttab BoboSS_hh_dc BoboSS_hh_oe using "$results\SummaryStatistics2.csv", replace cells("mean(fmt(3))" "sd(fmt(3) par)")


*Standard project: acreas=10, view=1, numboboter=2, NFledglings_recode=11
*******************************
*    SBC Model DHEU Table 5   *
*******************************
clear all
set more off
set seed 123456789
use "$dta\reanalysis_hurdle_final0.dta",clear


global X5 "acres view numboboter NFledglings_recode lnMinAmt PM UPA UPC age female mokids round2 Refund07 southto138 PObox lncoastdist"
global X9 "ppower age female dany denv moany momags moelect mokids round2 Refund07 southto138 lnparceldens PObox lncoastdist"
gen UPC_in_round2=UPC*round2
gen price_nonlinear=dcamt*(1-.5*exp(-.2*dcamt))/(.5-.5*exp(-.2*dcamt)) if dc==1
order dcamt price_nonlinear

 gen dcamt_sq=dcamt*dcamt

 drop if dcamt>offeramt&offeramt_all>0&dc==1
 replace offeramt=dcamt if dcamt<offeramt&offeramt_all>0&dc==1
 
 global varlist_ydc "acres view numboboter NFledglings_recode round2 Refund07 age female mis_indiv lncoastdist lnparceldens"
 global varlist_d "lnppower age female denv mokids mis_indiv lnparceldens lncoastdist round2"
 *global varlist_d "retcode lnppower age female denv mokids lncoastdist lnparceldens"

 global price "dcamt"
 sum $varlist_d
 drop if denv==.
 
 probit donate $varlist_d if dc==1
 mat list e(b)
 mat bparticipation=e(b)
 probit donate $varlist_ydc if dc==1&retcode==1
 mat list e(b)
 mat boutcome= e(b)
 mat list boutcome
 
 program define dhbobo_dc
    version 14
    args lnf beta1 beta2 beta6 beta3 beta4 beta5 
    tempvar d p non_l z p0 p1
    quietly gen double `d'=($ML_y1==1)
    quietly gen double `p'=normal(`beta1')
	quietly gen double `non_l'=(`beta3'-normal(`beta4'+`beta5')*`beta6')/(normal(`beta4'+`beta5')-normal(`beta4'))
	quietly gen double `z'=normal((`beta2'+`non_l'))
	*quietly gen double `z'=exp(`beta2'+`non_l')/(1+exp(`beta2'+`non_l'))  /*Logit*/
    quietly gen double `p0'=1-(`p'*`z')
    quietly gen double `p1'=(`p'*`z')
    quietly replace `lnf'=(1-`d')*ln(`p0')+`d'*ln(`p1')
 end 
 
 mat b=(bparticipation,boutcome,.01,.04,1,.144)
 ml model lf dhbobo_dc (hurdle: donate=$varlist_d)(above: donate=$varlist_ydc)(beta:dcamt,nocons)(warm_glow:)(mu:)(tao:dcamt,nocons) if dc==1,cluster(nameid)
 ml init b,copy
 ml maximize,iterate(10000) dif
 estat ic
 mat IC=r(S)
 eststo ModelDHEU
 estadd scalar AIC=IC[1,5]:ModelDHEU
 estadd scalar BIC=IC[1,6]:ModelDHEU
 esttab ModelDHEU using"$results\Model_DHEU.csv", replace se scalars(ll k AIC BIC) star(* .1 ** .05 *** .01) 
 
 scalar k_III=e(k)
 scalar LL_III=e(ll)
 
 
 predict Participant,equation(hurdle) 
 gen p_part=normprob(Participant)
 
 gen P=(p_part>=0.5)
 replace dcamt=0 if dc==0
 capture drop xb_donate alfa_diff C_star
 predict xb_donate,xb equation(above) 

 *Change to standard project level:  acreas=10, view=1, numboboter=2, NFledglings_recode=11
 gen alfa_diff=xb_donate+(_b[above:acres]*(10-acres)+_b[above:view]*(1-view)+_b[above:numboboter]*(2-numboboter)+_b[above:NFledglings_recode]*(11-NFledglings_recode))
 
 gen C_star=alfa_diff/(_b[beta:dcamt])
 gen WTP=C_star*p_part
 sum WTP if retcode==1
 sum WTP if donate==1
 
 capture drop nonl_cost_i main_i ll_III_i
 gen nonl_cost_i=(_b[warm_glow:_cons]-normprob(_b[mu:_cons]+_b[tao:dcamt]*dcamt)*_b[beta:dcamt]*dcamt)/(normprob(_b[mu:_cons]+_b[tao:dcamt]*dcamt)-normprob(_b[mu:_cons])) if dc==1
 gen main_i=exp(alfa_diff+nonl_cost_i)/(1+exp(alfa_diff+nonl_cost_i))
 gen ll_III_i=ln(1-p_part*(main_i)) if donate==0
 replace ll_III_i=ln(p_part*(main_i)) if donate==1
 egen LL_sum_i=sum(ll_III_i)
 
 save "$dta\reanalysis_hurdle_final1.dta",replace
 
use "$dta\reanalysis_hurdle_final1.dta",clear
*exponential 
gen W=_b[warm_glow:_cons]
gen tao=_b[tao:dcamt]
gen mu=_b[mu:_cons]
gen beta=_b[beta:dcamt]
gen WTO_cond=.
gen n=_n

*solving for WTO (willingness to donate/WTD in the paper): alfa_diff+(W-normprob(tao*WTO+mu)*beta*WTO)/(normprob(tao*WTO+mu)-normprob(mu))=0
*alfa_diff[1]+(W[1]-normal(tao[1]*a+mu[1])*beta[1]*a)/(normal(tao[1]*a+mu[1])-normal(mu[1]))
clear mata
set more off
putmata A=alfa_diff
putmata W=W
putmata tao=tao
putmata mu=mu
putmata beta=beta

mata:
void WTO( todo, p, a,w,t,m,b, v, g, H)
  {         
	        v = -(a+(w-normal(t*p+m)*b*p)/(normal(t*p+m)-normal(m)))^2
  }
Result=(1::5509)
for (i=1; i<=5509; i++) {
        a=A[i]
        w=W[i]
        t=tao[i]
        m=mu[i]
        b=beta[i]
        a

 S=optimize_init()

 optimize_init_evaluator(S, &WTO())
 optimize_init_evaluatortype(S, "d0")
 optimize_init_params(S, 1)
 optimize_init_argument(S,1,a)
 optimize_init_argument(S,2,w)
 optimize_init_argument(S,3,t)
 optimize_init_argument(S,4,m)
 optimize_init_argument(S,5,b)
 p=optimize(S)
 p
 Result[i]=p
} 

st_store(.,"WTO_cond",Result)
end
 
 gen WTO=WTO_cond*p_part
 order WTP WTO
 sum WTP WTO
 sum WTP WTO if P==1
 sum WTP WTO if retcode==1
 sum WTP WTO if donate==1
 
 hist WTP if donate==1
 hist WTO_cond if donate==1
 
 tab PM, sum(WTP)
 tab dc, sum(WTP)
  
 tab P dc,sum(WTP)
 twoway histogram WTP, bin(50) col(gs1) by(dc)

 tab P dc,sum(WTO)
 twoway histogram WTO, bin(50) col(gs1) by(dc)
 
 tab P dc,sum(WTO_cond)
 twoway histogram WTO_cond, bin(50) col(gs1) by(dc)
 *Simulation for revenue generated under DC

 set more off
 set seed 123456789
 capture drop simu_ER_* simu_Yes_*
 
 foreach p in 10 20 30 40 50 60 80 100 120 150 180 210 250{
 
 capture drop r y
 gen r=0
 gen y=0
 local N=2000
 forvalues i=1/`N'{
 capture drop simu_TR simu_TV u1 simu_P simu_Revenue simu_Donate
 gen u1=runiform()
 gen simu_P=(u1<p_part)
 gen simu_Donate=((WTO_cond*simu_P)>`p')
 gen simu_Revenue=`p'*simu_Donate
 egen simu_TR=sum(simu_Revenue)
 egen simu_TV=sum(simu_Donate)
 replace r=r+simu_TR
 replace y=y+simu_TV
 }
 gen simu_Yes_`p'=floor(y/`N')
 gen simu_ER_`p'=simu_Yes_`p'*`p'
 }
 capture drop simu_TR simu_TV u1 simu_P r y simu_Revenue simu_Donate
 sum simu_ER_*
 sum simu_Yes_*
 eststo Simu_DHEU_DC_R: estpost sum simu_ER_*
 eststo Simu_DHEU_DC_Yes: estpost sum simu_Yes_*
 esttab Simu_DHEU_DC_R Simu_DHEU_DC_Yes using "$results\Simu_DHEU_DC.csv",replace cells("mean(fmt(3))" "sd(fmt(3) par)")

 keep id nameid WTO WTP C_star WTO_cond
 save "$dta\predicted_WTO_WTP.dta",replace
 
 use "$dta\predicted_WTO_WTP.dta",clear
 sum WTO
 sum WTO_cond
 sum WTP 
 sum C_star

 
******************************
*   SBC Model DH - Table 5   *
******************************
set more off
use "$dta\reanalysis_hurdle_final1.dta",clear
global X5 "acres view numboboter NFledglings_recode lnMinAmt PM UPA UPC age female mokids round2 Refund07 southto138 PObox lncoastdist"
global X9 "ppower age female dany denv moany momags moelect mokids round2 Refund07 southto138 lnparceldens PObox lncoastdist"

 drop if dcamt>offeramt&offeramt_all>0&dc==1
 replace offeramt=dcamt if dcamt<offeramt&offeramt_all>0&dc==1

 global varlist_ydc "dcamt acres view numboboter NFledglings_recode round2 Refund07 age female mis_indiv lncoastdist lnparceldens"
 global varlist_d "lnppower age female denv mokids mis_indiv lnparceldens lncoastdist round2"
 
 sum $varlist_d
 drop if denv==.
 
 probit donate $varlist_d if dc==1
 mat list e(b)
 mat bparticipation=e(b)
 logit donate $varlist_ydc if dc==1
 mat list e(b)
 mat boutcome= e(b)
 mat list boutcome
 
 capture program drop dhbobo_dc_brief
 program define dhbobo_dc_brief
    version 14
    args lnf beta1 beta2 
    tempvar d p z p0 p1
    quietly gen double `d'=($ML_y1==1)
    quietly gen double `p'=normprob(`beta1')
    *quietly gen double `z'=exp(`beta2')/(1+exp(`beta2')) /*logit*/
	quietly gen double `z'=normprob(`beta2')    /*probit*/
    quietly gen double `p0'=1-(`p'*`z')
    quietly gen double `p1'=(`p'*`z')
    quietly replace `lnf'=(1-`d')*ln(`p0')+`d'*ln(`p1')
 end 
 
 mat b=(bparticipation,boutcome)
 ml model lf dhbobo_dc_brief (hurdle: donate=$varlist_d)(above: donate=$varlist_ydc) if dc==1,cluster(nameid)
 ml init b,copy
 ml maximize,dif
 
 estat ic
 mat IC=r(S)
 eststo ModelDH
 estadd scalar AIC=IC[1,5]:ModelDH
 estadd scalar BIC=IC[1,6]:ModelDH
 
 esttab ModelDH using"$results/Model_DH.csv", replace se scalar(ll k AIC BIC) star(* .1 ** .05 *** .01)
 
 scalar k_I=e(k)
 scalar LL_I=e(ll)
 scalar N=e(N)
 
 capture drop Participant p_part
 predict Participant,equation(hurdle) 
 gen p_part=normprob(Participant)
 
 capture drop P
 gen P=(Participant>0)
 replace dcamt=0 if dc==0
 capture drop xb_donate alfa_diff C_star
 predict xb_donate,xb equation(above)

 capture drop alfa_diff C_star 
 gen alfa_diff=xb_donate-_b[above:dcamt]*dcamt+(_b[above:acres]*(10-acres)+_b[above:view]*(1-view)+_b[above:numboboter]*(2-numboboter)+_b[above:NFledglings_recode]*(11-NFledglings_recode))
 
 gen C_star=alfa_diff/(-_b[above:dcamt])
 capture drop WTP
 gen WTP=C_star*p_part
 
 capture drop nonl_cost_i main_i
 capture drop main_i ll_I_i
 gen main_i=exp(xb_donate)/(1+exp(xb_donate)) if dc==1
 gen ll_I_i=ln(1-p_part*(main_i)) if donate==0
 replace ll_I_i=ln(p_part*(main_i)) if donate==1

 egen LL_sum_I_i=sum(ll_I_i)
 
 *Voung test
 scalar LR_star=LL_sum_i-LL_sum_I_i-log(N)*(k_III-k_I)/2
 di LR_star
 capture drop lli
 gen lli=ll_III_i-ll_I_i
 sum lli
 di r(Var)
 di r(N)
 scalar Vuong=LR_star/(sqrt(r(Var))*sqrt(r(N)))
 di Vuong
 di normal(Vuong)
 
 ****************************
 *   results in Table A1    *
 ****************************
 tab P donate
 tab P retcode
 tab P dc, sum(C_star)
 tab P dc, sum(WTP)
 twoway histogram WTP, bin(50) col(gs1) by(dc)
  
 set more off
 set seed 123456789
 capture drop simu_ER_* simu_Yes_*
 
 foreach p in 10 20 30 40 50 60 80 100 120 150 180 210 250{
 
 capture drop r y
 gen r=0
 gen y=0
 local N=2000
 forvalues i=1/`N'{
 capture drop simu_TR simu_TV u1 simu_P simu_Revenue simu_Donate
 gen u1=runiform()
 gen simu_P=(u1<p_part)
 gen simu_Donate=((C_star*simu_P)>`p')
 gen simu_Revenue=`p'*simu_Donate
 egen simu_TR=sum(simu_Revenue)
 egen simu_TV=sum(simu_Donate)
 replace r=r+simu_TR
 replace y=y+simu_TV
 }
 gen simu_Yes_`p'=floor(y/`N')
 gen simu_ER_`p'=simu_Yes_`p'*`p'
 }
 capture drop simu_TR simu_TV u1 simu_P r y simu_Revenue simu_Donate
 sum simu_ER_*
 sum simu_Yes_*
 eststo Simu_DH_DC_R: estpost sum simu_ER_*
 eststo Simu_DH_DC_Yes: estpost sum simu_Yes_*
 esttab Simu_DH_DC_R Simu_DH_DC_Yes using "$results\Simu_DH_DC.csv",replace cells("mean(fmt(3))" "sd(fmt(3) par)")

 keep id nameid WTP C_star
 save "$dta\predicted_WTP_dcdh.dta",replace
 *Note: this is actually WTD
 
 
 ************************************************
 * Results into Figures - Figure 2 and Figure 3 *
 ************************************************
 set more off
 use "$dta\predicted.dta",clear
 
 merge 1:1 id nameid using"$dta\reanalysis_bobo_hurdle.dta",keepusing(dc retcode retcode07)
 drop if _merge!=3
 drop _merge
 
 merge 1:1 id nameid using"$dta\predicted_WTO_WTP.dta"
 ren WTO WTO_ModelDHEU
 ren WTP WTP_ModelDHEU
 ren WTO_cond WTOcond_ModelDHEU
 ren C_star WTPcond_ModelDHEU
 drop _merge
 merge 1:1 id nameid using"$dta\predicted_WTP_dcdh.dta"
 ren WTP WTO_ModelDH
 ren C_star WTOcond_ModelDH
 drop _merge
 
 sum offer_hat if retcode==1
 local mean_offer_hat=round(r(mean),.001)
 kdensity offer_hat if retcode==1&offer_hat<200, xtitle("") title("Predicted OE offer (Mean=`mean_offer_hat')")  lw(thick) lc(gs0) sort
 graph save "$results\Model_OE.gph",replace

 
 sum WTO_ModelDH if retcode==1&dc==1
 local mean_WTD_DH=round(r(mean),.01)
 di `mean_WTD_DH'
 kdensity WTO_ModelDH if retcode==1&WTO_ModelDH<200, xtitle("") title("Model-DH DC WTD (Mean=`mean_WTD_DH')") lw(thick) lc(gs0) sort 
 graph save "$results\Model_DH.gph",replace

 sum WTO_ModelDHEU if retcode==1&dc==1
 local mean_WTD_DHEU=round(r(mean),.01)
 kdensity WTO_ModelDHEU if retcode==1&WTO_ModelDHEU<200, xtitle("") title("Model-DHEU DC WTD (Mean=`mean_WTD_DHEU')") lw(thick) lc(gs0) sort
 graph save "$results\Model_DHEU.gph",replace
 
 gr combine "$results\Model_DH.gph" "$results\Model_DHEU.gph", saving("$results\GraphII.gph",replace) imargin(tiny) cols(1) xsize(8) ysize(10) xcom ycom
 graph export "$results\\Figure2.tif",as(tif) replace
 
 capture drop offer_cond_g offer_hat_g WTOcond_ModelDHEU_g WTO_ModelDHEU_g WTPcond_ModelDHEU_g WTP_ModelDHEU_g WTOcond_ModelDH_g WTO_ModelDH_g 
 foreach V in offer_cond offer_hat WTOcond_ModelDHEU WTO_ModelDHEU WTPcond_ModelDHEU WTP_ModelDHEU WTOcond_ModelDH WTO_ModelDH {
sum `V'
local a=r(mean)

kdensity `V' if (`V'>0.001|`V'<-0.001)&(`V'>-500), gen(`V'_g) nograph at(`V')
}

twoway (connected WTO_ModelDH_g WTO_ModelDH if retcode==1, m(D) msize(vtiny) mcolor(gs5) lw(vvthin) lc(gs5) sort) ///
	   (connected WTO_ModelDHEU_g WTO_ModelDHEU if retcode==1, m(D) msize(vtiny) mcolor(gs0) lw(vvthin) lc(gs0) sort) , legend(order(1 "DC Model-DH" 2 "DC Model-DHEU"))

sum WTO_ModelDHEU, d
scalar WTO_p99=r(p99)
scalar WTO_p95=r(p95)
sum WTO_ModelDHEU if retcode==1&dc==1
local mean_WTD_DHEU=round(r(mean),.01)
count if WTO_ModelDHEU<0.01&WTO_ModelDHEU>-0.01& retcode==1&dc==1
local share_WTD_zero=round(r(N)/5509,.01)
di `share_WTD_zero'

sum WTP_ModelDHEU, d
scalar WTP_p99=r(p99)
scalar WTP_p95=r(p95)
sum WTP_ModelDHEU if retcode==1&dc==1
local mean_WTP_DHEU=round(r(mean),.1)
di `mean_WTP_DHEU'
count if WTP_ModelDHEU<0.01&WTP_ModelDHEU>-0.01&retcode==1&dc==1
local share_WTP_zero=round(r(N)/5509,.01)
di `share_WTP_zero'

twoway (connected WTO_ModelDHEU_g WTO_ModelDHEU if (WTO_ModelDHEU>0.01|WTO_ModelDHEU<-0.01)&WTO_ModelDHEU<100&retcode==1&dc==1, msize(vsmall) mcolor(gs0) lw(vvthin) sort) ///
	   (connected WTP_ModelDHEU_g WTP_ModelDHEU if WTP_ModelDHEU>-50&WTP_ModelDHEU<100&retcode==1&dc==1, msize(vsmall) mcolor(gs10) lw(vvthin) sort), saving("$results\GraphIII",replace) legend(order(1 "WTD(Mean=`mean_WTD_DHEU')" 2 "WTP(Mean=`mean_WTP_DHEU')"))
graph export "$results\\Figure3_1.tif",as(tif) replace

cap drop x s se
lpoly WTO_ModelDHEU WTP_ModelDHEU if WTP_ModelDHEU>-50&WTP_ModelDHEU<1000, gau nosc deg(2) bwidth(10) gen(x s) se(se) l(90) lineop(lw(thick) lc(gs0)) xtitle("WTP") ytitle("WTD") saving("$results\\Poly_WTO_WTP.gph",replace)
di r(bwidth)
graph export "$results\\Figure3_2.tif",as(tif) replace

